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Two defect particles that couple to a harmonic chain, acting as common reservoir, can become 
entangled even when the two defects do not directly interact and the harmonic chain is effectively 
a thermal reservoir for each individual defect. This dynamics is encountered for sufficiently low 
temperatures of the chain and depends on the initial state of the two oscillators. In particular, 
when each defect is prepared in a squeezed state, entanglement can be found at time scales at 
which the steady state of a single defect is reached. We provide a microscopic description of the 
coupled quantum dynamics of chain and defects. By means of numerical simulations, we explore 
the parameter regimes for which entanglement is found under the specific assumption that both 
particles couple to the same ion of the chain. This model provides the microscopic setting where 
bath-induced entanglement can be observed. 

PACS numbers: 03.67.Bg, 03.65.Yz, 05.40.Ca, 03.67.Mn 



I. INTRODUCTION 

It is commonly understood that the coupling of quan- 
tum systems to external environments destroys quan- 
tum effects, such as quantum superpositions and entan- 
glement. The microscopic picture is that this coupling 
generates correlations between the system and the en- 
vironmental degrees of freedom This results in an 
increase of the system's entropy, while the system state 
usually reaches a stationary state that is often well ap- 
proximated by a thermal state 0, Q . Such dynamics is 
well exemplified by the quantum Brownian motion, which 
can be microscopically modeled by the coupling of an os- 
cillator embedded in an ion crystal [^-Q- In Ref. [1], 
Rubin derived the conditions under which a defect oscil- 
lator thermalizes with the rest of the chain, which has 
been initially prepared in a thermal state at temperature 
T . This model provides an interesting realization of an 
Ohmic reservoir that contains in a natural way the rele- 
vant frequency scales. The physical system is closed and 
composed by one defect and the chain. From this per- 
spective it is important to mention most recent studies 
that analyze thermalization in closed systems j as 
well as recent proposals for simulating Ohmic reservoirs 
with chains of oscillators ■ 

Scaling up these dynamics by increasing the number 
of defects embedded in the crystal can lead to some sur- 
prises. Let us first assume that the parameters are chosen 
such that a single defect thermalizes with the rest of the 
chain. Contrary to the naive expectation that the two 
defects will reach a thermal state independent of their 
initial state, the two defects can be entangled by the 
reservoir at sufficiently low temperatures, even if they 
have been initially prepared in a separable state. This 
result can be ascribed to symmetries of the total Hamil- 
tonian that effectively decouple collective variables of the 
defect oscillators from the rest of the chain, leading to so- 



called decoherence free subspaces This mechanism 
of entanglement generation between two objects that are 
not coupled directly, but indirectly via a common larger 
physical system, has been discussed in various settings, 
see for instance [T2| - [25j . An important characteristic of 
most of these theoretical studies is the assumption that 
the two objects couple to an idealized bath with artifi- 
cially chosen spectral density. By contrast, in Ref. [2^ 
we considered the model of a one-dimensional harmonic 
crystal, whose spectral density was determined from ab 
initio calculations, and we showed that entanglement be- 
tween distant defects can be generated by the excitations 
of the crystal. However, a harmonic crystal does not 
always act as a perfect bosonic heat bath, since it can 
happen that the defect never relaxes to a steady state, 
even in the thermodynamic limit In the following, 
we perform a detailed investigation of the conditions un- 
der which (i) a generic harmonic chain plays the role of 
a thermalizing heat bath and (ii) two harmonic defects 
that couple to this chain are found to be entangled in the 
steady state. To this end, we numerically integrate the 
exact Heisenberg equations of motion of the total sys- 
tem, without making any weak-coupling or Markovian 
approximations. This allows us to explore the full pa- 
rameter regime. 

In this work we extend and complement parts of the 
findings reported in Ref. [1^ and systematically analyze 
the entanglement generation based on the microscopic 
model shown in Fig. [1] where the two defects couple at 
the same site of the ion chain. We examine the cor- 
relations between the defect oscillators for time scales 
that are smaller than the recurrence time (due to finite 
size effects), but for which a (quasi) stationary state is 
reached. Our objective is to connect our model predic- 
tions with previous studies on similar systems that were 
based on a phenomenological description of the reservoir 
[2l|-|2J] . For this purpose we tune the parameters 



2 




mi miN' 



FIG. 1: (Color online) Sketch of the microscopic model used 
for entanglement generation between the two defect oscilla- 
tors 1 and 2 (blue and red, respectively). The defects are 
confined by harmonic potentials with trap frequency Q, and 
couple with the strength 7 to a harmonic chain of A'' parti- 
cles. The particles of the chain interact via nearest-neighbor 
coupling with the strength k. The potentials with trap fre- 
quency ujb pin the edge oscillators of the harmonic chain. 



to a regime in which the chain effectively behaves like 
a (quasi) Ohmic reservoir. The numerical study allows 
us to determine both the stationary state, if it exists, as 
well as the out-of-cquilibrium dynamics for a vast range 
of parameters, for which a master equation description 
of the defect dynamics may not be convenient. The sim- 
ulations are supported by analytical investigations that 
yield a general criterion for the existence of steady-state 
entanglement. 

This paper is organized as follows: The microscopic 
model at the basis of our analysis is introduced in SecHIl 
Here, the basic idea leading to entanglement generation 
mediated by the chain is sketched. Section IIIII describes 
the theoretical formalism. The dynamics of the defects 
is studied in Sec. IIVI by means of a generalized quantum 
Langevin equation. The spectral density of the chain is 
discussed and the parameter regimes for which the har- 
monic chain acts as an Ohmic reservoir are identified. In 
Sec. |V] a detailed analysis of the entanglement behavior 
for different initial conditions and coupling parameters is 
given. The conclusions are drawn in Sec. lVIi and the Ap- 
pendixes |XHC] provide further aspects, as well as details 
of the calculations related to the Sec. IIIIIIVl 



strength 7 to the oscillator at one edge of the chain. The 
chain has been prepared in a thermal state at tempera- 
ture T . Our objective is to determine under which con- 
ditions the defect oscillators are entangled in the steady 
state. 



A. Hamiltonian 

The Hamiltonian determining the dynamics of the 
chain and the defect oscillators reads 



H = Hs + He 



Ht 



(1) 



and comprises the free Hamiltonian of the two defect os- 
cillators. 



p2 

2M 



2 



(2) 



the free Hamiltonian of the reservoir, 



N 



2m 



N-l 



i+lj 



and the interaction Hamiltonian, 



Hr 



(3) 



(4) 



wliich is assumed to be switched on at t = Q. 

Here, X^ denotes the position of the defect {n = 1,2), 
and Xi the displacement of the chain particle from the 
equilibrium position x^^^ = ia (i G {1, . . . A^}). With 
the corresponding canonically conjugate momenta 
and Pi, the nonvanishing commutation relations read 
[X^,P^] ~ \h and [xi,pi] = \h. Moreover, the slrorthand 
notation uji = usi^i^i + Si,]^) incorporates the trap fre- 
quencies of the edge oscillators in the chain. 



II. ENTANGLEMENT MEDIATED BY THE 
CHAIN 

In this section we first introduce the microscopic model 
that provides the basis of our study on entanglement gen- 
eration between two oscillators and then present the main 
idea why the two defect oscillators can become entangled 
via the interaction with the ion chain. 

The physical system is illustrated in Fig. [TJ It is com- 
posed of a chain of + 2 oscillators that couple with 
nearest-neighbor interaction. Among these, N oscillators 
have mass m and form an ordered linear chain with inter- 
particle distance a and interparticle coupling strength k. 
The oscillators at both ends of the chain are pinned by 
harmonic traps with frequency lub- The two additional 
defects have mass M and are confined by a harmonic po- 
tential with trap frequency f2. They couple with the same 



B. Basic idea of entanglement generation 

In presence of only one defect oscillator, the model in 
Fig. [H provides a generalization of Rubin's model [1]. Ru- 
bin showed in particular that the chain can act as a ther- 
mal bath for a single defect, provided some conditions 
are fulfilled, which involve the ratio M/m between the 
defect and the ions masses, the strength of the coupling, 
and the time scales in which the dynamics are analyzed. 
The scope of Sec. IIVI is to determine under which spe- 
cific conditions this dynamics is encountered for a finite 
chain. In this section we focus on the general idea and 
show that the ion chain can create entanglement between 
two defects, which are initially prepared in an uncorre- 
latcd quantum state. 

In general, bath-induced entanglement is endorsed by 
the symmetries of the Hamiltonian or, in the case of open 
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quantum systems, by the symmetries of the master equa- 
tion. We first observe that the total Hamihonian ([T|) is 
invariant under exchange of the coordinates of the two 
defect oscillators. It is therefore convenient to introduce 
center-of-mass (COM) and relative coordinates for the 
defect particles, 

X± = (Xi±X2)/V2, 

and the corresponding canonically conjugate momenta, 
P± = {Pi ± P2)/a/2, where the subscript + (-) de- 
notes the COM (relative) motion. In this representa- 
tion, the Hamiltonian ([l} can be written as the sum 
H = + H+, where 

governs the dynamics of the relative motion, and 




describes the coupling of the COM motion to the chain. 
Here, we denoted by 

= ^JWTlJM (7) 
the shifted trap frequency and by 

the chain Hamiltonian that includes the effect of the cou- 
pling constant 7 on the eigenspectrum. In this form it 
is evident that H- is a constant of motion: The rela- 
tive motion is decoupled from the chain. The COM, on 
the other hand, behaves as an effective defect particle 
that couples to one edge of the chain with the coupling 
strength •\/27. 

Under the conditions for which the chain acts as ther- 
mal bath for a single defect, it will induce thermalization 
of the COM defect particle and wash out possible ini- 
tial correlations between COM and relative motion of the 
defects. While the COM approaches a thermal state at 
temperature T after a transient time, the relative motion 
evolves freely and preserves some features of the initial 
states of the defects. 

The above-described dynamics is the key point in the 
creation of steady-state entanglement between the de- 
fects. For instance, if the relative motion is in a squeezed 
state and the temperature of the COM is sufficiently low, 
the product of the two orthogonal quadratures AX_AP+ 
(here taken in the reference frame rotating at the oscil- 
lator frequency il-y) can fall below the standard quan- 
tum limit giving rise to two-mode squeezing of the de- 
fects [13] and thus entanglement. The squeezing of the 
relative coordinate can be easily realized by preparing 
each individual defect in a squeezed state at the time 
t = 0. Figure m displays the contour plot of the logarith- 
mic negativity (28l . [29j that quantifies the entanglement 



between the defect oscillators. The logarithmic negativ- 
ity is shown as a function of the chain temperature T and 
of the initial squeezing parameter r of each defect oscilla- 
tor [22I, [l^l . The details of the calculations are provided 
in Sec. El 
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FIG. 2: Contour plots of the logarithmic negativity i5jv(f, T) 
as a function of the initial squeezing parameter f of the oscil- 
lators and of the temperature T = T/Ts, with Tg = fiflry/kB 
defining the temperature scale. The contour plot is evaluated 
for a chain of 1000 ions at times at which the COM motion 
has reached a stationary state, (b) Behavior at low tempera- 
tures. The black regions that are denoted by "NSD" indicate 
the parameter regime in which the defect oscillators are en- 
tangled in their steady state. The parameters are m = O.SAf, 
7 = 0.1k, and k = Mfi^. The squeezing parameters of both 
oscillators are taken to be equal. Further details are provided 
in Sec. |V] 

We note that this kind of dynamics has been predicted 
in Refs. [H,[13, where contour plots like the one in Fig. [5] 
have been introduced for the first time. In contrast to our 
work, the model used in Refs. [22,, 2M, takes advantage 
of the Hu-Paz- Zhang master equation [s^l and is based 
on a phenomenological treatment of the bath. In the 
present work, the bath is modeled by a chain of harmonic 
oscillators. Although we investigate a parameter regime 
in which our microscopic system reproduces the results 
of the Hu-Paz-Zhang master equation, we could likewise 
consider entanglement generation for a parameter regime 
in which the COM motion does not reach a thermal state. 
Such a regime, however, lies beyond the description based 
on the Hu-Paz-Zhang master equation [sol [3l|. 

We also would like to mention that entanglement me- 
diated by a chain of oscillators has been investigated in 
a series of works, such as [l5l l32| - [35| . In these works the 
chain is a homogeneous one-dimensional crystal and thus 
possesses discrete translational invariance. The regime is 
such that a unique stationary state exists in the thermo- 
dynamic limit which corresponds to a thermal state 
In Refs. [32| - [34| the authors characterize entanglement 
between two components of the chain in the steady state. 
The entanglement found in [l^ |35| between the ions at 
the chain edges is instead a dynamical effect, which ob- 
viously vanishes in the thermodynamic limit. 
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describes the thermal state of the chain at temperature T. 
Here, Z = Tr{exp {—(3Hb)} is the partition function and 
P = {ksT)^^ the inverse temperature with ks as Boltz- 
mann constant. Due to this choice of x(0)j there exist 
neither correlations between the defect oscillators nor be- 
tween the defects and chain at t ~{). 

More specifically, the defect oscillators are assumed 
to be prepared in pure states = \s^){s^\. Here, \s^) 
rewrite the reservoir Hamilto- denotes a squeezed state whose squeezing parameter 

s^^ = e"^^' is given by the absolute value > and 
the angle (pp, € (— tt, tt]. The corresponding first and sec- 



III. THEORETICAL FORMALISM 

In this section we develop the mathematical formal- 
ism, which we employ in the following sections for the 
characterization of the chain and for the analysis of the 
steady-state entanglement between the defects. 

For later convenience, we introduce the vector oper- 
ators for the reservoir particles x"^ = (xi, . . . ,2; at) and 

= (pi, . . . ,pn) and 
man ^ in the for m 

Hp 



with the potential matrix V G M given by 



P' 1 T 

2m 2 



(8) 



V = 



\ 



K — K 

K 2k 



2k 



\ 



+ kJ 



(9) 



The coupling between the oscillators and the reservoir 
induces a shift in the trap frequencies of the defect and 
chain particles, that depends on the coupling strength 
7. This effect can be highlighted by rewriting the full 
Hamiltonian ([1} in the form 



H = 



where 



2 

E 



p2 

2M 



M o o 



p_ 

2m 



ix^ 



7x1 {X1+X2) 



V + 2-iei®ei 



(10) 



(11) 



denotes the potential matrix including the shift due to 
the interaction. The quantity = (1, 0, . . . , 0) G is 
the first unit vector and (g) represents the dyadic product. 

An important point consists of the boundary condi- 
tions. For the model under consideration, we assume that 
the oscillators at both ends of the chain are confined by 
harmonic potentials with frequency lob- Although the 
potential of the ion at the other chain edge, j = N, has 
no influence on the dynamics of the defects for the time 
scales which are relevant to our analysis, we include it for 
symmetry reasons. As long as not specified elsewhere, we 
assume that lob = \/ n/m. throughout this paper. 



A. Initial states 



ond moments read {X^) = 

n 



(P.. 



and 



— 2r„ 2 , 2r„ • 2 

„ , , e cos -T" + e " sm 
2Mn V 2 2 



Km 



1 



2 iX,P, 
--sinh(2r^)sin0^ , 



(13) 



(14) 



(15) 



with (•) = Tr{-x(0)}. The moments in Eqs. p^-p^ 
define the initial covariance matrices (T^(0) of the defect 
oscillators at the time t = 0. 

For the following analysis it is also convenient to intro- 
duce the initial covariance matrix of the harmonic chain. 
We express it in terms of the individual block matrix 
elements 

CTxx(O) = (x ® x^) - (x) (g) {xf , 
app(0) = (p®p^)-(p)®(p)^ , 

fTxp(O) = i (x «) p^ + p «) x^) - (x) «) (p)"^ , 

whose explicit forms depend on the potential matrix Q 
and read m 



^xx(o) = -{mvy 



cothl'fr^ 



(16) 



^pp 



(0) = -(ml/) 2 coth 



"2" Im 



The initial state of the defect oscillators and the chain 
is given by the density matrix 

X(0) = pi ® P2 ® psiT) , 

where denotes the state of the defect oscillator 
{p = 1, 2) and 



Pb{T)^c^p{-I3Hb)/Z 



(12) 



together with o'xp(O) = and (x) = (p) = 0. 

B. Dimensionless variables 

With the total Hamiltonian and the initial covariance 
matrices at hand, we now introduce a dimensionless de- 
scription of our microscopic model. This reformulation is 
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system takes the rescaled form / 
satisfies the Heisenberg equation 



f{X^,P^;Xi,pi) and 



dl 
dt 



useful for the determination of the logarithmic negativity 
between the two defects. 

A typical length scale is the width of the ground state 
of the defect oscillator Hamiltonian ([2|), 

a-y = ^Jh/{^my) . 

The dimensionless position and momentum opera- 
tors for the two defects are defined as = X^/a^, 

Pf, = a^P^/h. For the oscillators^ of the reservoir we ac- matrix elements (HSJ read in' dimensionless form 
cordingly define Xi = Xi/a-y and Pi = a^pi/h. These def- 
initions imply the nonvanishing commutation relations 



We now come to the rescaled covariance matrices. 
With the dimensionless temperature T = ksT / [hVly), 
the inverse temperature /3 = ~ (3 hfl^, and the po- 
tential matrix V = V/{Mft'^), the nonvanishing block 



[Xt.Pi 



We further introduce the dimensionless mass m, trap 
frequencies Cjb and and coupling constants R and 7 
according to 

ffi — m/M ^ 

R = K/iMn"^) , 7 = 7/(Ml]2) 

With this choice, the mass of the defects M defines the 
unit mass, the shifted frequency il-y, Eq. ([7]), is the unit 
frequency, and the energy Mfi^ sets the relevant en- 
ergy scale. We note that the rescaled coupling strength 
7 = 7/(7 + MVi^) can only take on values in the inter- 
val < 7 < 1. Here, 7 = corresponds to 7 = 0, 
while 7 — > 1 represents the limit of infinitely large cou- 
pling 7 — > 00. 

The rescaled Hamiltonian H = H/{hil^) then reads 



^=1 



X. 



2 m 



2 



(17) 



7x1 



(^1 + X2 



with = V'''' /(Mfl^). The rescaled time is given by 
the variable 

For later convenience we also report the Hamiltonians 
governing the dynamics of relative and COM motion in 
their dimensionless form. They are given by 



H_ 



1 



(PI + Xt 



(18) 



lifhVyi cothf ^f- 
2^ ^ \2\m 



(20) 



1 / - -^ 1 

f^pp = 2^'^^)'' coth 



2 I TO 



Based on an appropriate one-to-one mapping r ~ r{r, (f) 
and (j) = 0(f, (j)) between the original and the new squeez- 
ing parameters > 0, 0^ £ (— 7r,7r], the covariance ma- 
trices for the defect oscillators p^ - p^ can be expressed 
in the convenient form 



1 



(21) 



In this expression, we introduced the symplectic and or- 
thogonal matrices (z G C) 



Siz) = 



z-^ 
z 



and 0(0) = 



(22) 

In this way, the elements of the initial covariance matrix 
for the defect oscillators (fT51)-(IT5]) reduce to 



. sin2 ^ 
2 



12 



■ sinh(2r^) sin < 



The above-mentioned one-to-one mapping is discussed in 
detail in Appendix [Xj The new parameters and (p^ de- 
fine the squeezing of the defect oscillators with respect to 
the shifted trap frequency fi^. Therefore, the squeezing 
parameter f = corresponds to the ground state of a 
harmonic oscillator with trap frequency VL^. 



and 



H+ = - [Pi + xi 



p 

2fh 



F<^' X - X+ (7^x) , 
(19) 

where we have introduced the dimensionless coupling vec- 
tor 7"^ = (\/2 7,0,...,0) e R^. 

According to these definitions, an operator function 
f{Xfj_,Pfj_;Xi,Pi) acting on the Hilbert space of the total 



C. Formal solution of the equations of motion 

The formal solution of the Heisenberg equations of mo- 
tion for the position and momentum operators of both 
defect and bath oscillators can be written as a linear 
map between their initial and final values. For this 
purpose, we introduce the vector of the position and 
momentum operators for defect and chain oscillators. 
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= (Xi,Pi,X2,P2,x^,p^_) e E,4+2^, and rewrite the 
total Hamiltonian P7)) as H = ^ C^Ti, C, with the posi- 
tive definite matrix H. Furthermore, we introduce the 
antisymmetric block matrix 



J 



J2 
Jn 



-J' 



-J- 



that contains the submatrices 



10 0' 
-10 
1 
0-10, 



and Jjss = 



In 

-Ijv 



(23) 

Here, In & M,^'^^'^ denotes the identity matrix. 

By means of these definitions, the Heisenberg equa- 
tions of motion for the position and momentum operators 
reduce to 



dt 



Their formal solution reads 
with the symplectic matrix 



(24) 



m 



„JHt 



The time evolution of the total covariance matrix, V(t), 
is given in terms of the linear mapping by the relation 



v(t) = r(t)v(o)r^(t). 



(25) 



V(0) = 



where V(0) is the total covariance matrix at t — 0, 
which is composed of the initial covariance matrices (pO)) 
and fill) and takes the form 



/(Tl(0) 

(72(0) 


V 



Equation ((25|) represents the basis of the numerical sim- 
ulations used in the analysis of entanglement generation. 
In this context, the covariance matrix of the defect oscil- 
lators 5](F) is of particular interest. It is extracted from 
the total covariance matrix V{t) according to 




[S(i)],, = [V(t)]. 



(26) 



with i,k € {1, 2, 3, 4}. Since we aim at the determination 
of the steady-state entanglement, it suffices to evaluate 
the covariance matrix I](i) at times t > tth- Here, tth 
represents the time scale at which the COM defect oscil- 
lator reaches a stationary state, provided the harmonic 
chain acts as a thermal bath. For this reason, we examine 
in the next section the conditions for which the reservoir 
displays this behavior. 



IV. CHARACTERIZATION OF THE 
RESERVOIR 

The harmonic chain plays a basic role in our study of 
entanglement generation between the defects for the fol- 
lowing reason: Although the total dynamics is unitary 
and the system is finite, the chain can act as a thermal 
bath for the COM motion of the defects, while the rela- 
tive motion is uncoupled. In order to understand under 
which conditions this mechanism leads to entanglement, 
a detailed knowledge about the action of the chain on the 
COM motion is necessary. Hence, the purpose of this sec- 
tion is to characterize the chain in terms of a reservoir 
and identify the parameter regime for which these condi- 
tions are valid. 



A. Generalized Quantum Langevin Equations for 
the defects 

Let us consider the dynamics of the defect oscillators. 
The dynamics of the relative motion is governed by the 
Hamiltonian ([T5)) . and the solution of the corresponding 
Heisenberg equations of motion simply describes the evo- 
lution of a harmonic oscillator with frequency that 
reads 



X^{t) = X^{0)cost + P_ (0) sin t , 
P-{t) = -X_(0)sini-f P_(0)cosi, 



(27) 



where we recall that t = Qjt. The COM motion, nev- 
ertheless, remains coupled to the oscillator at the chain 
edge. We rewrite its equation of motion in terms of a 
generalized quantum Langevin equation (GQLE). Start- 
ing from the Heisenberg equations of motion for the op- 
erators X+, P_|_, X, and p, the GQLE follows by formal 
integration of the equations for the chain degrees of free- 
dom and takes the form 



dp 



V(t-0 



dX^ 

~dF 

(t)X+(0) 



dt' + (i-T+{o))x+{t) 



(28) 



Here, we have introduced the memory-friction kernel [5[, 
which reads 



N 



(29) 



for t > 0, while it vanishes otherwise. We have also 
introduced the operator-valued random force, which is 
defined by [sj 



N 



7+cos(w+f)x+(0) 



3in(cS+t)p+(0) 



(30) 
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In the expressions for the memory-friction kernel and the the corresponding frequency obtained in the hmit 



random force, the quantities uij' and 7^ appear. The first 
ones denote the eigenfrequencies of the potential matrix 
V<-^'> given by Eq. ([TT|) . They follow from the diagonal- 
ization of the chain potential V'-''' and are defined by the 
relation 



ol v^-'^o+ 



m ■ diag((ajj^ 



+ \2 



('^^f), (31) 



where 0+ is the orthogonal matrix which diagonalizcs 
V^''K In particular, the orthogonal matrix 0_|_ estab- 
lishes the relation between the normal and the ori gina l 
coordinates, x+ = O^x and p+ = O^p, see e.g. j3a |. 
The quantities x+(0) and p+(0) in Eq. ^ stand for 
the j-th component of the vectors x+ and p+, respec- 
tively. The parameters 7^ are the coupling strengths 
to the j-th normal mode of the reservoir and are given 
by 7_|_ = 7. In the following we adopt the conven- 
tion that the eigenfrequencies are ordered according to 
< < < . . . < u;^. 

An important quantity that characterizes the influence 
of the reservoir on the COM motion is the environmen- 
tal spectral density. This quantity is the Fourier cosine- 
transform of the memory-friction kernel (|29p 



J+(w) 



cj I r+(t) cos(wF) dt 
'0 

- ^ it'' 



(32) 



The spectral density provides important insight into the 
action of the chain on the dynamics of the COM motion 
of the defect. 

Before we proceed, we characterize the chain's nor- 
mal modes. The eigenfrequencies are the solutions of 
Eq. (|5T|) . which includes the shift due to the coupling 
of the defects with the edge ion. By appropriately set- 
ting the frequency of the edge potentials to the value 
lob = y^fi/m (see Appendix IB|) . we obtain for the nor- 
mal mode spectrum in the limit 7 — > 



sm 



kjtt 



(33) 



where to 



(0) 



UJ^'^^kj) and kj ^ jTr/{a{N + 1)) is the 



wave number, which appropriately denotes the modes 
when the Bloch theorem applies and takes on the val- 
ues (j G {1,...,A^}). This expression agrees with the 
one found for periodic boundary conditions The fre- 
quency (Dcut — y^4:K/m is the high-frequency cutoff. The 
resulting spectrum, Eq. (p3)) . is displayed in Fig. |3l 
The eigenmodes, however, are, strictly speaking, not 
phononic waves. 

Let us now consider the case 7 > 0. We expect for suffi- 
ciently small 7 that the effect of this coupling on the chain 
normal-mode spectrum is negligible. To quantify this 
statement, we consider the difference ACuj = uij' — 
that involves the eigenfrequency cD^ given by Eq. and 



7—^0. Figure [3ljb) displays the corrections Awj for dif- 
ferent coupling strengths 7 = 0.1,0.15,0.2 and constant 
K = 1. For these values, the difference Alj^ is much 
smaller than all other physical parameters. 




0.5 
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FIG. 3: (Color online) (a) Spectrum ajj"' of the potential V, 
Eq. @, for the parameter values m = 0.5, 7 — 0.1, and R = 1 
and (b) corrections Aujj to the eigenfrequencies of 1/'^' for 
the coupling strengths 7 — 0.1 (solid), 7 — 0.15 (dotted), and 
7 = 0.2 (dashed line). The chain is composed of A*' = 1000 



Figures Hlja) and IH^b) display the spectral density for 
a choice of the parameters 7 and R, and taking fh = 0.5. 
For most of the parameter values the spectral density is 
linear about the value w = 1, which corresponds to the 
frequency of the defect oscillator. In this case, the chain 
acts like a (quasi) Ohmic environment. A change in the 
mass ratio ifi affects the spectral density in so far as the 
eigenfrequencies scales with ~ leading to a 

change in the bandwidth Hicnt- 



(a) ^ 



(b) 

10'V+(w) 



K = 1.0 



7 = 0.1 



FIG. 4: (Color online) Spectral density J+(ijj) as a function 
of Cj. The parameters in (a) are R = 1, while 7 — 0.1 (solid), 
7 = 0.15 (dotted), and 7 — 0.2 (dashed line). In (b) we 
take 7 = 0.1 while n = 0.5 (dashed), k = 1 (solid), and 
K = 1.5 (dotted). In both cases the mass ratio is m = 0.5. 
The spectral densities is linear in the vicinity of the oscillator 
frequency to = 1. A slight difference is found for the cases 
7 = 0.2, K = 1 (a) and 7 = 0.1, R = 0.5 (b) where J+{uj) ~ ui° 
with a > 1 in the vicinity of tj = 1. 



B. Thermodynamic limit 

In Ref. [3] Rubin showed that a chain of oscillators 
with one embedded defect, exhibiting a spectrum as in 
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Eq. (|33|) . can act like an Ohmic bath for the defect par- 
ticle. This behavior is found in Rubin's model provided 
that the temperature of the chain is finite, the mass ra- 
tio satisfies m/M < 1, and the thermodynamic limit is 
taken, which corresponds to the limit N ^ oo keeping 
the interparticle distance a in the chain constant [37l | . Fi- 
nite size effects are found for times larger than the time 
scale trcvj which is discussed in the next section. How- 
ever, they can be neglected if (i) the defect oscillator 
reaches a stationary state over time scales ith such that 
tth ^ trev and (ii) the analysis can be restricted to these 
time scales. 

We now discuss these assumptions in relation to our 
microscopic model, where, different from Rubin's model, 
the coupling strength 7 appears in addition to the cou- 
pling constant k. We are specifically interested in identi- 
fying the parameter regimes for which the effective defect 
of our model, the COM, thermalizes with the rest of the 
chain. 

For this purpose we first consider the formal solu- 
tion of the Heisenberg equations of motion in Eq. (p4)) 
in terms of the linear mapping 7+(f) =e'^+'^+* (here 
the positive definite matrix and the antisymmetric 
matrix J7+ are defined in analogy to the discussion of 
Sec. IIII C|) . The GQLE can formally be solved by apply- 
ing a Laplace transformation to both sides of Eq. (pS]) . 
which yields an algebraic equation for the Laplace trans- 
form of X^(t). In this case the residue theorem can be 
applied ;38| : The simple poles of the integrand are deter- 
mined from the eigenfrequencies of the positive definite 

~ ~l/2~ ~l/2 

matrix W+ = T|_ V+ T|_ with the block matrices 

r+=(ji) and ^+=(4 7<^')' ^^^^ 

which respectively originate from the kinetic and poten- 
tial energy part of the Hamiltonian in Eq. ([6]). The 
sum of the residues yields a quasiperiodic function which 
is equivalent to the expression for Xj^(t) deduced from 
the linear mapping T+{t). As in Rubin's model of a single 
defect in a one-dimensional crystal thermalization is 
found in the limit iV — >■ 00 due to the formation of a con- 
tinuous frequency band, provided no isolated frequencies 
above the frequency band occur. The existence of such 
isolated frequencies would result in residual oscillations 
of the coupled defect (in our case the COM motion) at 
long times. 

Here, we arc interested in the parameter regime in 
which the COM motion of the defect reaches a stationary 
state before the time scale frcv, and in particular, when 
this stationary state is a thermal state at the tempera- 
ture T in which the chain was initially prepared. In order 
to identify the coupling strengths 7 and k for which this 
is verified, we perform a numerical search of the isolated 
frequencies of the matrix W+ The results are presented 
in Figs. [S]for mass ratios to = 0.5 and to = 1. In the re- 
gion above the broad curves no isolated frequency of W+ 
was found. In this domain thermalization of the defect 



COM occurs in the thermodynamic limit according to 
our numerical simulations. In the segmented region be- 
low the boundary curve, at least one isolated frequency 
of W+ exists. Here, the labels on the contour lines in- 
dicate the value of the largest isolated frequency of the 
normal modes. The dots in the upper region indicate the 
parameter values used in our simulations: They all lie in 
the region where no isolated frequencies exist. 



(a) (b) 




7 7 



FIG. 5: (Color online) Diagram illustrating the existence of 
isolated frequencies of the matrix W+ for different coupling 
strengths 7 and R and the two mass ratios fh = 0.5 (a) and 
m = 1 (b). At least one isolated frequency is found for the 
parameters below the black boundary line: The value of the 
largest one is indicated by the contour lines. The dots in the 
white region indicate the parameter values used in our simula- 
tions: They all lie in the region where no isolated frequencies 
exist. 



C. Finite chains 

Since our analysis is essentially numerical, we consider 
finite chains and we aim at observing a (transient) sta- 
tionary state of the COM motion before finite size effects 
become relevant. The latter are characterized by the time 
scale ficv = L/cs, where L = Na is the chain length and 
Cg — (Iicut a/2 the sound velocity [11]. The time scale ?rev 
grows linearly with the particle number N, showing that 
by choosing a sufficiently large number of particles ther- 
malization of the defect particle could be observed before 
finite-size effects become significant (which we denote by 
"revivals"). 

We illustrate the thermalization of the COM mo- 
tion of the defects by showing the time evolution 
of the variances AX^(t) = (X^(t)) - (X+(f))^ and 

APlit) = (Plit)) - {P+{i)y in Fig. H After a tran- 
sient time, the variances approach their stationary val- 
ues AXl{f) and AP^(T) that depend on the initial tem- 
perature of the chain, but not on the initial squeezing 
of the defects. The appearance of revivals after ij-cv is 
clearly visible. 

A good estimate for AXl{T) and APl(T) fol- 
lows from the assumption that the total system (de- 
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FIG. 6: Time evolution of the variances AX^(f) (black) 
and AP+(f) (grey) for the temperature T — 10~^ of the chain. 
After a transient time, the covariance matrices reach the sta- 
tionary values AX^(f) = 0.5031 and AP^(f) = 0.4988, 
which indicate that the harmonic chain causes squeezing of 
the momentum variable at very low temperatures. The dot- 
ted vertical line represents the revival time fth ~ 1416. The 
other parameters were chosen to be m = 0.5, k = 1, 7 = 0.1, 
<i>\=<i>2 = 0, and fi = f2 = i ln(l - 7). 



fects and harmonic chain) is in the thermal state 
Pth = exp(— /3/f)/Zth. Here, H denotes the total Hamil- 
tonian (|17[) . /? the inverse temperature of the initial 
chain (|12p. and Zth ~ Tr{exp(— /3ff)} the corresponding 
partition function. With the help of Eqs. ([20]) and (p4)) . 
we find the following values for the stationary state vari- 
ances of the defects: 



coth 



APlif) = - 



wi 



coth I 



The indices on the right-hand side of the equations indi- 
cate the (1, l)-elements of the matrices inside the brack- 
ets. The fact that this estimate works so well, despite 
the unitary time evolution of the total system, is remi- 
niscent of the concept of "canonical typicality" 0, Q that 
recently gained a lot of attention. 

The parameter values of our microscopic model have 
to meet several constraints. First of all, the coupling 
strength R must be sufficiently large in order to guar- 
antee that the frequency of the two oscillators lies well 
below the cutoff frequency Wcut and more specifically in 
the linear region of the spectral density. Moreover, the 
value 7 must be sufficiently small such that the disper- 
sion spectrum of the harmonic chain is not significantly 
perturbed by the coupling with the defects. There is also 
a further bound to the coupling strength 7 that stems 
from the necessity to reduce computational resources. In 
fact, 7 determines the rate at which the center-of-mass 
motion reaches a stationary state. Very small values of 
7 would require that one chooses an increasing particle 
number N in order to observe a stationary state well 
before fiov, which results in a formidable computational 
problem. 

In order to account for all these requirements, we have 



chosen the parameter values m = 0.5, 7 = 0.1, and 
K = 1 as standard parameters for our numerical simula- 
tions throughout this paper. As in the previous figures, 
we illustrate the changes in the numerical results that 
arise from different coupling constants by using the two 
parameter sets: (i) The ^-variation parameters. The re- 
sults are presented for three different coupling strengths 
7 = 0.1, 0.15, 0.2 and for the fixed parameter value R = 1. 
(a) K-variation parameters. The results are presented for 
constant 7 = 0.1 and for variable k = 0.5,1,1.5. For the 
case 7 = O.l, K = 1.5 we used N = 2000 ions in the 
chain, while in all other cases it was sufficient to work 
with N = 1000 ions in order to observe that the COM 
motion reached a stationary state well before the revival 
time ircv 

It is instructive to analyze the variances of the COM 
position and momentum, after the steady state has been 
reached. Figure [7] shows the variances AX^(T') and 
AP^(f) for a large temperature range T S [0,6] (a) 
and for small temperatures T <S [0, 0.2] (b) given the 
"/-variation parameters with rh = 0.5. For large tem- 
peratures (a) the variances grow linearly with a slightly 
different slope for the individual coupling strengths 7, 
whereby AX^{f) > AP^{T). In the low temperature 
regime (b), we observe squeezing of the COM momen- 
tum AP^(T') < 1/2 that increases for larger coupling 
strengths 7. We note that this squeezing is induced by 
the coupling with the bath and has been identified in 
the studies reported in [2^ [2^ . It is reminiscent of the 
squeezing found for large coupling strengths in the Drude 
model [alii. 



(a) 

AXl{T),APl(,T) 




FIG. 7: (Color online) Variances AX^(T) (three upper 
curves) and AP^{T) (three lower curves) of the COM mo- 
tion after thermalization for large temperatures T £ [0, 6] (a) 
and for low temperatures T G [0, 0.2] (b). The parameters are 
m = 0.5, K = 1, and 7 — 0.1 (solid), 7 = 0.15 (dotted), and 
7 = 0.2 (dashed line). 

Figure [5] displays the corresponding behavior of the 
COM variances for the R-variation parameters. As be- 
fore, we find a linear behavior of AX^ (T) and APl{f) 
for large temperatures (a). For low temperatures (b), the 
squeezing of the variance AP^(T') < 1/2 becomes larger 
for smaller R and vice versa. 
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(b) 

AXl(T),APl{T) 
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FIG. 8: (Color online) Variances AX^(T) (three upper 
curves) and /S.P'^{T) (three lower curves) for large and low 
temperatures, T € [0,6] (a) and T £ [0,0.2] (b), respectively. 
The parameters are m = 0.5, 7 = 0.1, and R = 0.5 (dashed), 
R = 1 (solid), and R — 1.5 (dotted). 
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FIG. 9: (Color online) Memory-friction kernel r+(t). The 
plots in (a) and (b) are evaluated for the parameters of the 
curves for the spectral density in Fig. [4] (a) and (b), respec- 
tively. 



V. CHAIN-MEDIATED ENTANGLEMENT 
BETWEEN THE DEFECTS 



Finally, we emphasize that a mass ratio m = 1 leads 
only to marginal changes in the temperature behavior of 
the variances AX^ (T) and APl{f). 

These properties directly affect the behavior of bath- 
mediated entanglement between the defect oscillators, as 
we show in Sec. |Vl 



D. Memory effects 



In this section we report the predictions of our model 
on correlations between the two defect oscillators. Entan- 
glement is quantified by means of the logarithmic nega- 
tivity [28I [23 , that is evaluated using the covariance ma- 
trix at the time scale at which the COM motion of the 
two defects has reached a (quasi) steady state (before the 
revival time). We present the results for the logarithmic 
negativity found for different choices of parameters, such 
as the initial squeezing of the defect oscillator, the tem- 
perature of the chain, and the coupling strength between 
chain and defects. 



We now analyze memory effects in our model using 
our parameter choice. For this purpose we discuss the 
memory- friction kernel ^+{t) of the GQLE that is 
connected to the spectral density J+(w) for < > accord- 
ing to relation 



f+(t) 



2 

I" ./o 



which inverts Eq. (|32p . For strict Ohmic dissipation, 
the memory- friction kernel would read r_|_(<) — 2T 5{t) 
with r as friction constant, and the GQLE would reduce 
to the ordinary Langevin equation without memory ef- 
fects, provided that the "slip term" ~2T 5(t) Xj^[Q) and 
the oscillator frequency shift —2T 5{0) X+{t) can be dis- 
regarded [Hj- However, in our model we do not meet 
the requirements of a strict Ohmic environment since the 
cutoff frequency Wcut is only a few times larger than the 
oscillator frequency of the two coupled oscillators. For 
this reason, memory effects are present. 

The figures [9ja) and (b) display the memory-friction 
kernel as a function of time: An oscillatory decay is ob- 
served over a time scale t that is of the order of one, 
corresponding to t ~ l/fi-,, [4l|. Hence, non-Markovian 
effects arc present, but irrelevant for the dynamics of en- 
tanglement generation between the defects, as is shown 
in the following. 



A. Logarithmic negativity of the oscillators 

Since the state of the two defects is a Gaussian state 
at all times, the most convenient entanglement measure 
for our purpose is the logarithmic negativity [H, . 

In what follows, we present the final results and refer to 
Appendix [U] for further details on the calculations. The 
logarithmic negativity reads 



-Ejv(tth) = max{0,£:jv(ith)} 



(35) 



and contains the function £'Ar(fth) = — ln(2 i'_(fth))j 
which depends on the smallest symplectic eigen- 
value u-{tt\i) of the partial transpose of the covariance 
matrix I](ith), defined in Eq. (|26p. The covariance ma- 
trix E(fth) describes the state of the system for suffi- 
ciently long times, f ~ fth, after which the COM motion 
has reached a stationary state independent of its initial 
state. The smallest symplectic eigenvalue v-{tt\^) follows 
from the identity [3 SI 



th) 



A2(ith)-4detE(ith) 



(36) 
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with the timc-indcpcndcnt determinant 
1 



det S(tth) = z 1 + cosh(2fi) cosh(2f2) 

8 V 

- cos(A0) sinh(2ri) sinh(2f2)) , (37) 

and the osciUating auxiUary function 

A(tth) = Ao + A2 cos(2tth + (^) . (38) 

In the last two expressions, we have introduced the rela- 
tive squeezing angle A0 = 4>2 ~ <t>i^ well as the coefB- 
cients 

Ao = i(AX2 +Ap2^(cosh(2fi)+cosh(2f2)) (39) 



and 



A2 



AXl - AP2 



sinh^(2fi) + sinh^(2f2) 
2 cos(A0) sinh(2fi) sinh(2f2))' . (40) 



The constant phase Cp can be determined, but is of no fur- 
ther interest to us. Due to the periodicity of the auxiliary 
function (|38p the quantity £N{ith) oscillates for < i-^^v 
between a minimal and maximal value and £'^^^. 

The formulas provide a generalization of pre- 

viously obtained expression for the logarithmic negativ- 
ityHm. 

Following the nomenclature of [22, 1231 j distinguish 
three qualitatively different situations for the entangle- 
ment of the two oscillators, (i) When < 0, the 
logarithmic negativity is zero and we find no entangle- 
ment between the oscillators. We call this scenario the 
sudden death (SD) phase because any transient entangle- 
ment disappears abruptly before the thermalized state 
is reached, (ii) When < < f we obtain an 
alternating sequence of periods with zero and nonzero 
logarithmic negativity, the so-called sudden death and 
revival (SDR) phase, (iii) Finally, when > the 
two oscillators are entangled after thermalization which 
we call the no sudden death (NSD) phase. In Fig. [TU] 
we exemplify these different phases by showing the time 
evolution of £N{tth) for three initial temperatures. Fig- 
ure [TUlJa) displays the long-time behavior of f7v(ith), its 
evolution toward the steady state. Here, the occurrence 
of revivals after trov are visible. Figure llOf b) zooms in 
the behavior at ith ~ O.Otrcv, showing that the loga- 
rithmic negativity exhibits oscillations at the frequency 
fL,- These oscillations have been also observed in Rcfs. 
[23 . [23I and their physical origin simply lies in the de- 
coupling of the relative coordinate from the rest of the 
dynamics. In fact, the squeezed variance of the relative 
motion rotates with frequency 20.^, and correspondingly 
the smallest symplectic eigenvalue oscillates at the same 
frequency. We also note that, by choosing the squeezing 
parameters according to fi = r2 = 0, we find by virtue of 
Eqs. dMl) and (gOl) that the logarithmic negativity ([35]) of 



the steady state becomes time independent and displays 
no further oscillations. The underlying reason is that for 
fi = f 2 = the initial state of the relative motion corre- 
sponds to the ground state of the Hamiltonian ([5]) . 
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FIG. 10: Time evolution of the logarithmic negativity £N{i) 
for three different temperatures T = 10~^, 0.27, 0.33 (from top 
to bottom). The curves exemplify the different behaviors of 
the entanglement (NSD, SDR, SD). The dotted vertical line in 
(a) represents the revival time frov. (b) Detail of the behavior 
about the time i « O.Ofrov The other parameters are m = 
0.5, K = 1, 7 = 0.1, cj>i = (f)2 — and f\ = = \ ln(l — 7). 
Here, frov « 1416. 



In this context, we would like to point out that the 
NSD phase can be characterized by the fulfillment of the 
inequality (see also Eq. ([C5|) in Appendix [C|) 



Ao - A2 - 4 det S(fth) - ^ > , 



(41) 



which follows from £7v(ith) = — lii(2 i'-(ith)) > or 
equivalently i>_(fth) < \ evaluated for the minimal value 



of A(tth), Eq. ([381). Thus, if inequality ([IJ) is satis- 
fied, the two defect oscillators are entangled after the 
COM has reached a stationary state (before the revival 
time frev)- This inequality in connection with the iden- 
tities (|57|) - (|M|) provides a general criterion for the exis- 
tence of steady-state entanglement for arbitrary initial 
squeezed states of the defects. 



B. Entanglement generation for different initial 
parameters and coupling strengths 

In this section we report the logarithmic negativity of 
the defect oscillators after the COM defect oscillator has 
reached a stationary state, for different values of the ini- 
tial squeezing of the defects and of the initial temper- 
ature of the ion chain. The results are displayed using 
the type of contour plots first introduced in Ref . [23| , 
which highlight the different entanglement regions (NSD, 
SD and SDR) as a function of the modulus of the initial 
squeezing parameter and the temperature of the reser- 
voir. 

We first consider the case in which the initial states 
of the defect oscillators are characterized by the same 
squeezing parameters, f\ = f2 = f and A(/) = 0. We 
use the inequality AX^(T) > APl{f), which we found 
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numerically for the considered parameter regime, and re- 
duce the determinant ([57]) to the form 



detE(tth) - -AXlAPl, 
while the coefficients (P^l) and PO]) read 



Ao =i (^AX^ + APl^ cosh(2f) 
A2 =i (ax2 - AP|) sinh(2f) . 



(42) 

(43) 
(44) 



These expressions lead to the following simple form of the 
entanglement condition (|4T|) that characterizes the NSD 
phase: 



.2 2f 



AP+e 



+ AX2 , 



-2r 



AX^ AP; 



1 

4>0- 



With the substitution y = e^'', this relation reduces to 
a quadratic inequality in y that yields two independent 
conditions for the steady-state entanglement of the two 
oscillators. The first of these conditions reads 



(45) 



and tells us that entanglement between the oscillators 
will occur at any temperature T as long as the initial 
squeezing parameter f is sufficiently large. The underly- 
ing mechanism for this entanglement generation is based 
on the existence of a decoherence free subspace, following 
from the decoupling of the relative motion. 

The second entanglement condition takes the form 



APl{f)<^e- 



(46) 



and is only satisfied for sufficiently small squeezing pa- 
rameters f and temperatures T. We call this second 
mechanism bath-induced entanglement because it arises 
from the squeezing of AP^ (T) caused by the interaction 
of the oscillators with the reservoir. It is clear that the 
two mechanisms are competing. 

Figure [5] displays the different phases of entanglement 
for varying r and T including the contour lines of the 
logarithmic negativity within the NSD region. In Fig. 

one can observe the behavior at large squeezing and 
high temperatures. Here, entanglement in the NSD re- 
gion is due to the decoupling of the relative motion and 
is determined by the condition (P5|). The SDR region 
is not visible, but lies between the NSD and SD phases. 
Figure [2Jb) displays the behavior at small squeezing and 
low temperatures. One can here observe the NSD island, 
which occurs in the vicinity of f = 0, T = and is sep- 
arated by the SDR phase from the main NSD region. 
This island stems from the bath-induced entanglement 
according to Eq. (|46)) . 

Since the squeezing of the COM motion at low tem- 
peratures is rather small, the NSD region due to bath- 
induced entanglement covers only a small region of 



Fig. [U^b). The size of the region can be increased by in- 
creasing the squeezing of the variance AP^{T). Accord- 
ing to Fig.[71[b), this can be achieved by increasing the pa- 
rameter 7. Figure [TT] displays the corresponding contour 
plots in the regime of small squeezing parameters and low 
temperatures for two values of the coupling strength 7: 
An increase of the NSD region of bath-induced entan- 
glement is observed for larger coupling strengths 7. We 
recall, however, that this behavior can saturate, when 7 
takes values at which the transient steady state is not 
reached before trcv The squeezing of the COM vari- 
ance can also be increased by decreasing the coupling 
strength R, as illustrated in Fig. IHJb). Figure [T2] depicts 
the change in the entanglement behavior for varying R. 
Here one can see that the size of the region where bath- 
induced entanglement is found is larger for smaller values 
of K. 
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FIG. 11: Contour plots of the logarithmic negativity i5jv(r, T) 
for m = 0.5, k = 1, and (a) 7 = 0.15, (b) 7 = 0.2. 
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FIG. 12: Contour plots of EN{f,T) for m = 0.5, 7 = 0.1, and 
(a) K = 0.5, (b) R = 1.5. 



When the two oscillators are instead prepared in 
squeezed states with a relative squeezing angle A0 7^ 0,, 
the entanglement will be diminished. In fact, such initial 
states lead to a smaller squeezing of the relative motion. 
A representative situation is found for A0 = tt, namely, 
when the squeezed quadratures of the defect oscillators 
are orthogonal. In this case the relative motion is not 
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squeezed and one obtains for the determinant pT]) and 
the coefficients p9p and ([Ml with fi=f2=f the ex- 
pressions 

detli(ith) = ^ AX^ AP| cosh^(2f) , 

as weU as 

Ao = ^ (^AXl + APl^ cosh(2f) and A2 = . 

In this case, the entanglement condition (|^T1) reduces to 

^(^AXl+APl^ cosh{2f)-AXlAPl cosh^(2f)-i > 0. 

This inequahty is fully equivalent to the new criterion 

1 



AP_f (T) < 



2 cosh(2f) ' 



which is only satisfied for a squeezed COM momentum, 
in analogy to Eq. p6|) . It shows that entanglement be- 
tween the defects can only be generated by the active 
coupling with the bath. The existence of a decoherence 
free subspace does not support entanglement generation 
in this case. 

Thus, the relative squeezing angle A(j) can be used as a 
control parameter to distinguish between the two mech- 
anisms that lead to steady-state entanglement. This ob- 
servation makes our model a favorable microscopic set- 
ting to study the generation of bath-induced entangle- 
ment. 



VI. CONCLUSIONS 

We have numerically investigated the dynamics of two 
defects coupled to one edge of a harmonic crystal and 
identified the parameter regime for which the defect vari- 
ables reach a quasi steady state. This (quasi) equilibrium 
sets in for time scales which are smaller than the re- 
vival time scale characterized by finite size effects. From 
its features and its scaling behavior for different system 
sizes, we can conclude that it corresponds to the equi- 
librium reached in the thermodynamic limit, when the 
number of ions of the chain is infinitely large. The anal- 
ysis of the correlations between the defects shows that 
they can become entangled in the steady state. Such en- 
tanglement emerges as a consequence of the symmetries 
of the Hamiltonian, and it follows the dynamics outlined 
in Refs. j22l . [23| where it was determined by means of 
an effective master equation mimicking the effect of the 
bath. Our work shows that a physical model, such as 
the considered extension of Rubin's model, establishes 
a microscopic realization of this dynamics. It allows us 
to determine the relevant time scales which emerge from 
the spectral properties of the chain, the defects, and their 
mutual coupling. Moreover, it gives us the possibility to 
analyze the dynamics in regimes where a master equation 



approach is not convenient (e. g. when finite size effects 
become relevant). 

This work provides a microscopic understanding of the 
dynamics of bath induced entanglement, building upon 
the general criterion given by Eqs. ([35 |) -(|4T |) . Based on a 
realistic model, it goes beyond the reach of idealized set- 
tings studied so far that employ ideal bosonic heat baths 
with artificially chosen spectral densities. An interest- 
ing next step will be the extension of our model to non- 
Gaussian initial states and nonquadratic Hamiltonians 
for the defects. As long as the symmetry is preserved, we 
anticipate that the underlying mechanisms will support 
the formation of steady-state entanglement. Whether 
such an extension will lead to an enhancement in the 
entanglement generated between the defects is however 
an open question. 

In a follow up to this article we will discuss the gen- 
eration of entanglement between two defects that couple 
to distant sites of the chain, thereby extending and com- 
plementing the findings reported in Ref . [2^ , which were 
not addressed in the present article. 
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Appendix A: Transformation of the squeezing 
parameters 



Based on the dimensionless description of Subsec llII Bl 
and the matrices S{z) and 0{(j)), Eq. ((22)) . we find for the 
initial covariance matrices of the defects cr^(O), Eqs. p^ - 
P^. the dimensionless form 

(7^(0) = i S0)O^{^^) S{e^^^) 0{cj>^) S0) . (Al) 

Due to the outer symplectic matrices <S'(f2), we would 
arrive at much more complicated expressions for the 
logarithmic negativity in Sec. IV Al when starting from 
Eq. (|A1[) . These expressions would conceal the class of 
squeezing parameters that lead to the same steady-state 
entanglement between the defects. 

For this reason, it is advantageous to introduce new 
squeezing parameters = r^e"^'' that overcome these 
difficulties by transforming the covariance matrices (jAl[) 
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to the simpler form (j21[) . The corresponding transforma- 
tion equations r = r(f, (p) and = (f)[f, (f>) follow directly 
from the diagonalization of ()A1|) and a subsequent com- 
parison of the resulting eigenvalues and eigenvectors with 

Eq. dm). 

In this way, we find that the one-to-one mapping be- 
tween the new squeezing parameters > 0, (/>/i £ (— tt, tt] 
and the original ones (r^ > 0, 0^ € (— tt, tt]) depends 
on n and splits into three different domains of defini- 
tion. Since we have < 1, the mapping r — r{f, (j)) and 
4> = 4>{f, 4>) reads for the special case = (f > 0) 



r(f,0) 
0(f,O) 



(A2a) 



For = TT (r > 0) we find accordingly 



r(r, tt) 



ilnrj) • sign(r-h ilnrj) , (A2b) 



0(f,7r) = TT - e (f + ilnfl) , 

where Q{x) denotes the Heaviside step function. 
The mapping of the remaining open domain f > 0, 
(j) e (-TT, tt) \ {0} onto r > 0, e (-tt, tt) \ {0} is finally 
given by the one-parameter family of coordinate trans- 
formations 



r(r, 0) = - arcosh {R+) 



(A2c) 



0(r, 4>) = 2 arctan 



IR\-1-R^+ sinh(2f) sin0 



iR\-l + R-+ sinh(2f) sine 



with the auxiliary functions R± = R±(r, 0, 51) defined by 



1 -\ sinh(2f) 



^^^\cosh(2f) 

n 2 



Substitution of the transformation Eqs. (jA2[) into the 
original covariance matrix (jAl[) yields directly the con- 
venient form ((2T|) . 

The inverse transformation equations f = r(r, 0) and 
= 0(r, 0) follow in analogy to Eqs. (|A2[) by simply 
replacing the role of < 1 in the derivation with its 
inverse f2~^ > 1. Again, the domain of definition splits 
into three different parts. For the special case = 
(r > 0) the inverse mapping is given by 



fir, 0) = (r + i In n) ■ sign (r + i In il) 



(r,0) = 7r(l-ei 



For = TT (r > 0) it reads accordingly 



r(r, n) = r - 
4>{r, tt) ~ TT . 



■In 17, 



(A3a) 



(A3b) 



As above, we find for the mapping of the open domain 
r > 0, e (-TT, tt) \ {0} onto f > 0, € (-tt, tt) \ {0} a 



slightly more complicated expression 



1 



r(r, 0) = - arcosh 



(A3c) 



-2 arctan 



'R^-l-R^- sinh(2r) sine 



iR^-l + R-- sinh(2r) sine 



with the auxiliary functions R± = R± (r, 0, ^ ) given by 



- 1 \ sinh(2r) 

Si =F ^ ^ — - cos ( 

n 2 



Using these inverse transformations, one can determine 
the values of the new squeezing parameters f ^ and 0^ for 
a given set of initial squeezing parameters and 0^. 

The effect of the transformation Eqs. (jA2[) is illustrated 
in Fig. [12] by showing the coordinate lines r = r{f, 0) 
and = 0(f, 0) for fj = i and constant values of f 
and 0. The rather small value for Cl was only chosen to 
highlight the effect of the transformation Eqs. (jA2[) . and 
does not correspond to any of the parameter values used 
throughout the paper. 




FIG. 13: (Color online) Illustration of the curvilinear coordi- 
nates defined by Eqs. (|A2|I for H = i. The red curves indicate 
the coordinate lines r = r{f,<j)) and = (p{f,<j)) for constant 
f > and varying G (— vr, tt) \ {0}. The blue curves are 
obtained for constant and varying f. 

We conclude this appendix by pointing out that the 
transformation Eqs. (jA2p arc only valid for < 1. In 
the case of f2 > 1, one obtains the corresponding trans- 
formation equations by simply interchanging Eqs. (|A2[) 
and its inverse ()A3|) and replacing r o f and o 0. 



Appendix B: Spectral density for different trap 
frequencies Cjb 

The purpose of this appendix is to show that the shape 
of the spectral density depends crucially on the choice of 
the edge frequency Qb- In the main part of the paper, we 
restrict ourselves to the fixed value ujb ~ y^K/fh. In this 



15 



way, we compensate for the missing frequency shift of 
the ions at the end of the chain (they couple only to one 
neighboring ion). This choice yields a suitable tridiago- 
nal form for the potential matrix ^ whose eigenvalues 
and eigenvectors can be anal ytic ally determined using 
the methods outlined in Ref. [4al- As a result, we find 
the spectrum in Eq. p3p for the specific trap frequency 

Since the two defect oscillators couple to the edge par- 
ticle of the harmonic chain, the trap frequency ujb has 
an immediate infiuence on the behavior of the reservoir. 
In order to illustrate this fact, we show in Fig. [HT a) the 
spectral density for the standard parameters 7 = 0.1, 
K = 1, m = 0.5 and the three different trap frequen- 
cies ujb = 0.1 (dashed curve), uib = (solid curve), 
and = 2 (dotted curve). Whereas J+{oj) exhibits a 
pronounced non-Ohmic behavior for small trap frequen- 
cies ljb ^ V2, it still displays a linear growth in the 
neighborhood of w = 1 for > y/2. 



(a) 

10^ J+(w) 



(b)_ 

io'r+{f) 








8 



12 



t 



FIG. 14: Spectral density (a) and the memory- friction ker- 
nel (b) for the parameters 7 = 0.1, R, — 1, fh = 0.5 and 
the three different trap frequencies lDb = 0.1 (dashed line), 
UJB ~ (solid line), and ljb = 2 (dotted line). 



tivity for arbitrary initial squeezing parameters and 
steady-state variances AX^ (T) and APl{T). 



Logarithmic negativity and covariance matrix in 
COM and relative coordinates 



We start by recalling the definition of the combined 
vector of the position and momentum operators for the 
two defect oscillators ^ = (Xi, A; ^2, P2V ■ The corre- 
sponding covariance matrix is given by the expression 
Sa^j = i (^ + e/3_Ca) - (Ca) i^p) with a, /3 e {1, . . . 4}. 
It can be rewritten in the block form, 



A C 
B 



(CI) 



where A^B £ E,^, denote the covariance matrices of the 
first and second defects, respectively, and C G E,^ char- 
acterizes the correlations between them. Next, we define 
the partially transposed covariance matrix S = A S A 
with the help of the diagonal matrix A = diag(l, 1,1,-1). 
The logarithmic negativity [1^, [2^ can then be deter- 
mined from the smallest symplectic eigenvalue v- of S 
and reads 



En = max{0, -ln(2!>_)} . 



(C2) 



We note that the symplectic eigenvalues of E coin- 
cide with the common, positive eigenvalues of the ma- 
trix — iJ2S, where J2 is given by Eq. (|23p . 

It is possible to write down an explicit expression for 
the smallest symplectic eigenvalue [3, and for this 
purpose, we introduce the function 



Figure I14r b) depicts the corresponding memory- 
friction kernel r+(f). For Qb = 0.1 we find a nonoscil- 
latory, slowly decaying function r+(t) which indicates 



large memory effects in the GQLE (f28| . For ujb ~ 2 
we obtain an oscillatory behavior of the memory-friction 
kernel; however, the oscillations do not decay for large 
times. The reason for this behavior is the existence of an 
isolated frequency in the spectrum of W+ which prevents 
the COM motion from thermalization, see Sec. IIVI 



A(S) = det A det B 2 det C 

which is invariant under symplectic transformations. By 
applying this function to the partially transposed covari- 
ance matrix S, we obtain the auxiliary function 

A A(AEA) = dctA + detS - 2detC (C3) 

With this quantity at hand, the smallest symplectic 
eigenvalue of S follows from the identity 



Appendix C: Analytic expressions for the 
logarithmic negativity of the steady state 

In this appendix we derive the analytic expressions 
used in Sec. |V] for the evaluation of the logarithmic neg- 
ativity. We first recall how to find the logarithmic nega- 
tivity for a given covariance matrix in general 'Q^2^M,. 
We then rewrite this formalism in COM and relative co- 
ordinates and apply it to the specific covariance matrix of 
the defects after they reached the steady state. Finally, 
we sketch the derivation of the simple expressions p7p . 
([M)) and PO)) that provide the logarithmic nega- 



VA2 -4detE 



(C4) 



Given the covariance matrix in block form (|C1[) . we thus 
determine the logarithmic negativity (|C2p by evaluating 
the smallest symplectic eigenvalue ()C4p with the help of 
the auxiliary function (|C3[) . 

Entanglement between the two oscillators is only found 
when Ej<s = — \w{2vJ) > which is equivalent to !>_ < i. 
Using Eq. (jC4[) . one can show that this entanglement 
condition coincides with the Simon criterion 14711 



A-4detI]- T > 0. 

4 



(C5) 
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Now, due to the decoupling of the relative coordi- 
nate of the two defect oscillators in our microscopic 
model, we seek for an expression of £>_ that is based 
on the covariance matrix in COM and relative coor- 
dinates. For this reason, we define in analogy to 
above the combined vector for the COM and rela- 
tive coordinates = P+; P_)^ and write 
^ap = f(e<*'C^*; for the cor- 

responding covariance matrix with block form 



(±) 



(C(±))T 5(±) 



(C6) 



With the transformation matrix 



R 



1 

71 



I2 I2 
I2 -I2 



R-y 



the connection between the COM and relative coordi- 
nates and their corresponding covariance matrices reads 



|<±'=i?| and S< 



(C7) 



In order to rewrite the quantities that appear in the 
smallest symplcctic eigenvalue (|C4[) in terms of the block 
matrices A^^\ 5'=*=' and (7'="=', we take advantage of the 
fact that the transformation matrix R is symplectic. An 
immediate consequence of this observation is the validity 
of the identities 



det S = dot S 



(±) 



(C8) 



and 



A(S) = A(i? R^) = A(E'±') . (C9) 
With the help of Eq. (|C9p , we easily find the relation 

A = A(E<±')-det[A<±'-B<±' + (C<±'f-C'±'] . (CIO) 



In conclusion, the smallest symplectic eigenvalue (jC4[) , as 
well as the entanglement condition (jCSp can be directly 
determined from the covariance matrix in COM and rel- 
ative coordinates (jC6|) by means of the identities (|C8[) 
and (ICTOj) . 



2. The covariance matrix after thermalization of 
the COM motion 

The manifestation of correlations between the defects 
is a direct consequence of the decoupling of the relative 
coordinates and the thermalization of the COM motion. 
This statement can be well illustrated my means of the 
covariance matrix of the defect oscillators. Initially, the 
covariance matrix of the two defects reads 



S(0) = 



<Tl(0) 
^2(0) 



where the 5'^j(0) are given by Eq. ([2T|) . The transforma- 
tion to COM and relative coordinates via Eq. (|C7p yields 
the covariance matrix 



I]<^'(0) = - 



1 f ai(0) + CT2(0) ai(0)-a2(0) 



Si(0)-CT2(0) ai(0)-fa2(0) 



which displays correlations between the COM and rela- 
tive coordinates as long as the initial squeezing parame- 
ters of the two defect oscillators differ. 

After turning on the coupling to the reservoir, the 
COM motion of the two defects thcrmalizes after a tran- 
sient time tth < frev which gives rise to the covariance 
matrix 



(7+(T) 
a- (tth) 



(CU) 



Here, the time-indcpcndcnt submatrix of the COM reads 
a+(T)=('^5 A*?.^') (C12) 



and contains the variances 



API 



AXl 



AXliT) 



and 



AP?(T) on the diagonal. The actual values of 



AA"^ and AP^ are numerically determined and depend 
on the initial temperature T of the reservoir. The time- 
dependent covariance matrix of the relative coordinate 



-(tth) =T_(tth)(T-(0)TZ^(ith) 



(C13) 



describes the free time evolution of the initial covariance 
matrix 



^"(0) = ^(ai(0)+a2(0)) 
with the help of the orthogonal matrix 



T_(t-) = 



cos t sin t 
- sin t cos t 



(C14) 



(C15) 



that follows from the solution (P?]). By transforming 
the covariance matrix (jClip back to the original coor- 
dinates, we finally obtain the covariance matrix of the 
steady state. 



S(t 



_ 1 / a+{f) + a-{tth) (T+(r)-a-(tth) 
- 2 I a+{T) - a- (tth) a+{T) + a" (tth) 



which now exhibits correlations between the first and sec- 
ond defect oscillator. 



3. Derivation of the auxiliary functions for the 
logarithmic negativity 

In this section, we present the main steps of the deriva- 
tion of the analytic expressions (|57)) and ([55)1 . which are 
used for the evaluation of the logarithmic negativity in 
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Sec. IV Al We thereby take advantage of the determinant 
identity 



det{A ± B) = det A + det B 
± dot 



(C16) 



( an 


ai2 


)±detr" 


&12 \ 


[ bn 


b22 , 


/ V "21 


022 / 



which holds true for any two matrices A = (a.ifc) G C^^^ 
and B = (bik) e C2x2. 

In order to find the expression for the determinant p7l) , 
we first recall Eqs. (pT]) . (p^ and ((CT3)) to ob- 

tain 



det S(ith) = AX^ API det 5-^(0) , 



(C17) 



Using the definition of (7^(0), Eq. ()C14p . and inserting 
the initial covariance matrices of the two oscillators (piT) , 
we find with Acj) = 02 — 01 

det(T-(0) = Y^det [^(e2^i) + O^(A0)5(e2^^)O(A0)] . 

The last expression can be easily evaluated with the help 
of identity (|C16|) . which yields after some minor algebra 

deta"(0) =i(l + cosh(2fi)cosh(2f2) 
8 V 

- cos (A0) sinh(2fi) sinh(2f2)) . (C18) 



Substitution of the last expression into Eq. (|C17|) pro- 
vides the expression (|37p for the determinant of S(tth)- 
To derive the time-dependent auxiliary function psp . 
we start from Eqs. (jClOp . (jClip . and (|C13p and obtain 
with the orthogonality of r_ (t) 

A(ith) =deta+(ith;T) 4-detCT-(0) 

-det[a+(tth;T)-a-(0)] , (C19) 



where 



By applying the identity (|C16|) to Eq. (|C19p , we find after 
a straightforward calculation 



A(tth) - -{AXl + API) (ari(O) + ^22(0)) 

+ ^{AXl ~ API) [(ari(O) - ^22(0)) cos(2fth) 
+ 2c7r2(0)sin(2fth)] . 

When we combine the two terms in the bracket of the 
last equation in a single cosine, we obtain the general 
form ([55)) of the auxiliary function. The resulting co- 
efficients Aq and A2 can be rewritten in terms of the 
determinant and trace of a~ (0) according to 



Ao 



and 



AXl + APl)Tr{a- (0)} 



(C20) 



A2 = . 



AXl ~ API 



V(Tr{a-(0)})2 -4detCT-(0). 

(C21) 

From the initial covariance matrices (|2ip and the defini- 
tion (|C14[) . we obtain for the trace 



Tr{CT"(0)} = -( cosh(2fi) -Hcosh(2f2; 



which together with the determinant (jClSP finally yields 

(Tr{5--(0)})2-4deta-(0) = i ( sinh^(2fi) + sinh2(2f2) 
-h2cos(A0) sinh(2fi)sinh(2f2)^ 



a+(ith; T) = T^itti,) a+if) T_(tth) • 



Substitution of the last two expressions into Eqs. (jC20[) 
and (jC21[) finally concludes our derivation of the coeffi- 
cients ^ and (gni). 
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